# -*- cording:utf-8 -*-

# scipy & numpy & plot Libraries
from scipy import *
from pylab import *
from enthought.mayavi import mlab
# from Standard Library

# from User Defined Library
import util 

# ======================================================
# class section
# ======================================================
class Index(object):
    """prepar like arrays of structure for index Matrix
    """
    class P:
        pass
    class Ux:
        class plus:
            pass
        class minus:
            pass
    class Uy:
        class plus:
            pass
        class minus:
            pass
    class Uz:
        class plus:
            pass
        class minus:
            pass
    def plot(self):
        fig = figure()
        ax = fig.add_subplot(111)
        ax.hold(True)
        ax.plot(index.P[:,0], index.P[:,1], "s",  label="P.field")
        ax.plot(index.Ux.plus[:,0]+.5, index.Ux.plus[:,1], "r>", label="Ux.plus")
        ax.plot(index.Ux.minus[:,0]+.5, index.Ux.minus[:,1], "r<", label="Ux.minus")
        ax.plot(index.Uy.plus[:,0], index.Uy.plus[:,1]+.5, "g^", label="Uy.plus")
        ax.plot(index.Uy.minus[:,0], index.Uy.minus[:,1]+.5, "gv", label="Uy.minus")
        ax.plot(index.Ux.field[:,0]+.5, index.Ux.field[:,1], "r.", label="Ux.field")
        ax.plot(index.Uy.field[:,0], index.Uy.field[:,1]+.5, "g.", label="Uy.field")
        ax.legend()
        ax.grid(True)
        show()

class Num(object):
    """prepar like arrays of structure for index Matrix
    """
    def __init__(self):
        self.X = 0
        self.Y = 0
        self.Z = 0

class Room(object):
    def __init__(self, sizeMat):
        sizeMat = list(sizeMat)
        self.X = float(sizeMat[0])
        self.Y = float(sizeMat[1])
        if len(sizeMat) is 3:
            self.Z = float(sizeMat[2])
        else:
            self.Z = None
        
class Field(object):    
    def __init__(self, _size, _dxyz, _temperature=15.0, _rho=1.21, _bound="wall"):
        """
        roomX:Room Size X-direct.
        roomY:Room size Y-direct.
        fs: sampling Frequency
        temperature: temperature of air
        """
        util.attributesFromDict(locals())
        self._c = 331. + 0.6*_temperature # set soundSpeed
        self._kapper = _rho * self._c**2
        print self._kapper
        # Set Boundary Condition of Matrices
        # Very very easy setting. Truth be told, for cleary use, I'll have to set as reactance.
        if _bound is "air":
            self._matrixBoundCondtion = _c*_rho #
        elif _bound is "wall":
            self._matrixBoundCondtion = 10000 
        
        self.Z = 1000
        self.setFieldFromdxyz()
                
    def setFieldFromdxyz(self):
        num = Num()
        if self._size.Z is None:
            num.X = int(self._size.X/self._dxyz[0]) + 1
            num.Y = int(self._size.Y/self._dxyz[1]) + 1
            self._dxyz = [self._size.X/num.X, self._size.Y/num.Y]
            # last "/1.1" makes "dt" safety. >> ref:Courant Condition (CFL Condition)
            self._dt = 1./(self._c*sqrt((1./self._dxyz[0])**2 + (1./self._dxyz[1])**2))/1.1
            self._fs = 1./self._dt
            if self._c*self._dt <= 1./sqrt(2./self._dxyz[0]**2):
                print "courant condition ok!!"
                print self._dt/(self._dxyz[0]*self._rho)
            else:
                print self._dt
                print dxyz
                raise ValueError, "check !!"
        else:
            num.X = int(self._size.X/self._dxyz[0]) + 1
            num.Y = int(self._size.Y/self._dxyz[1]) + 1
            num.Z = int(self._size.Z/self._dxyz[2]) + 1
            self._dxyz = [self._size.X/num.X, self._size.Y/num.Y, self._size.Z/num.Z]
            # last "/1.1" makes "dt" safety. >> ref:Courant Condition (CFL Condition)
            self._dt = 1./(self._c*sqrt(self._dxyz[0]**2 + self._dxyz[1]**2 + self._dxyz[2]**2))/1.2
            self._fs = 1./self._dt
        self._num = num
        msize = [num.X+2, num.Y+2]
        self.P1 = zeros(msize, "float64")
        self.P2 = zeros(msize, "float64")
        self.Ux1 = zeros(msize, "float64")
        self.Ux2 = zeros(msize, "float64")
        self.Uy1 = zeros(msize, "float64")
        self.Uy2 = zeros(msize, "float64")

class IRRecorder(object):
    def __init__(self, Field, IRList):
        self._Pnow = Field.P2
        self._fs = Field._fs
        self._list = IRList

    def recIR(self):
        for ir in self._list:
            ir.rec(self._Pnow)
    
    def saveIR(self, dirName, bits=32):
        ii = 1
        for ir in self._list:
            if is_string_like(dirName):
                util.chkMakeDir(dirName)
                filename = "./" + dirName + "/IR%3d.wav" % ii
            else:
                raise ValueError, "Check director name."
            ir.outputWav(self._fs, bits, filename)
            ii += 1

class IR(object):
    def __init__(self, FDTD, position):
        self._position = FDTD.pos2ind(position)
        self._record = []
        
    def rec(self, pressure):
        self._record.append(pressure[self._position[0], self._position[1]])
      
    def outputWav(self, fs, bits, filename):
        util.wavwrite(self._record, fs, bits, filename)

class FDTD(Field):
    def setSource(self, type, sourcePos):
        self._sourcePos = self.pos2ind(sourcePos)
        if type == "gaus":
            tau0 = 48./self._c*min(self._dxyz)
            tau = arange(0, 3*tau0, self._dt)
            alp = (4./tau)**2
            self._source = exp(-alp*(tau-tau0)**2)
            #plot(self._source, "o--"); show()
        elif type == "sine":
            pass
        else: 
            raise NameError, "Not correct source type"
            
    def chkCalcField(self, nodes, index):
        """check the calculation area
        matrix"formX" (X is P/Ux/Uy) has 0/1. "0" means calc, "1" means NOT calc indices there. 
        """
        print " check calculation field start. ".rjust(50)
        num = self._num
        matSize = [num.X+2, num.Y+2]
        points = array([[x,y] for x in linspace(0,self._size.X,num.X) for y in linspace(0,self._size.Y,num.Y)])
        inIndex, onIndex, outIndex = util.inpoly(points, nodes)
        inIndex = array(list(set(r_[inIndex, onIndex])))
        
        # Prepare logical(0/1) matrix. "1/-1" means calculation points. 
        tmp = zeros(points.shape[0])
        tmp[outIndex] = 1
        tmp = tmp.reshape([num.X, num.Y])
        P = zeros(matSize)
        P[1:-1,1:-1] = tmp
        UxObj = vstack([diff(-P,axis=0), zeros([1, P.shape[1]])])
        UyObj = hstack([diff(-P,axis=1), zeros([P.shape[0], 1])])
        UxField = P * (1-abs(UxObj))
        UyField = P * (1-abs(UyObj))

        index.P = util.ind2sub(matSize, find(P==1))
        index.Ux.plus = util.ind2sub(matSize, find(UxObj==1))
        index.Ux.minus = util.ind2sub(matSize, find(UxObj==-1))
        index.Ux.field = util.ind2sub(matSize, find(UxField==1))
        index.Uy.plus = util.ind2sub(matSize, find(UyObj==1))
        index.Uy.minus = util.ind2sub(matSize, find(UyObj==-1))
        index.Uy.field = util.ind2sub(matSize, find(UyField==1))
        print " end. ".rjust(50)


    def looper(self, index, recorder, calcTime):
        self.chkDim()
        print " loop calculation start.... ".rjust(50, "=")
        s = mlab.surf(self.Ux2[1:-1,1:-1])
        cb = mlab.colorbar()
        ttl = mlab.title("sec: %4d[ms]" % 0)
        print "source calculation start.... ".rjust(50)
        
        for cnt in range(self.chkCalcTimes(calcTime)):
            if cnt < len(self._source):
                print "%3d::: source:%.3f,  P1:%.2f".rjust(60) % (cnt, self._source[cnt], self.P1[self._sourcePos[0], self._sourcePos[1]])                
                self.P1[self._sourcePos[0], self._sourcePos[1]] += self._source[cnt]
                print "afterP1:%.2f".rjust(60) % self.P1[self._sourcePos[0], self._sourcePos[1]]                
            else:
                print "souece calculation end.... ".rjust(50)            
            self.updater(index)
            self.plotter(s, cb, ttl, cnt)
            recorder.recIR()
        print "loop calculation end.... ".rjust(50)
        
    def chkDim(self):
        exec("self.updater = self.update" + str(len(self._dxyz))) 

    def chkCalcTimes(self, calcTime):
        return int(calcTime/self._dt)

    def update2(self, index):
        # update P =====================================
        for ind in index.P:
            ii = ind[0] 
            jj = ind[1] 
            if (ind==self._sourcePos).all():
                print "sourcePos P2 :%.4f, P1:%.4f".ljust(20) % (self.P2[ii,jj], self.P1[ii,jj])
            self.P2[ii,jj] = self.P1[ii,jj] \
                             - self._kapper*self._dt*((self.Ux1[ii,jj]-self.Ux1[ii-1,jj])/self._dxyz[0]\
                                                      +(self.Uy1[ii,jj]-self.Uy1[ii,jj-1])/self._dxyz[1])
            if (ind==self._sourcePos).all():
                print "sourcePos P2A:%.4f, P1:%.4f".ljust(20) % (self.P2[ii,jj], self.P1[ii,jj])
        # update Ux  ===============================
        for ind in index.Ux.field:
            ii = ind[0]
            jj = ind[1]
            if (ind==self._sourcePos).all():
                print "sourcePos U2 :%.4f, U1:%.4f, check:%f".ljust(20) % (self.Ux2[ii,jj], self.Ux1[ii,jj], self._dt/(self._dxyz[0]*self._rho) )
            self.Ux2[ii,jj] = self.Ux1[ii,jj] - (self._dt/(self._dxyz[0]*self._rho)) * (self.P2[ii+1,jj] - self.P2[ii,jj])
            if (ind==self._sourcePos).all():
                print "sourcePos U2A:%.4f, U1:%.4f".ljust(20) % (self.Ux2[ii,jj], self.Ux1[ii,jj])
                #print self.Ux2
            
        for ind in index.Ux.plus:
            ii = ind[0] 
            jj = ind[1]
            self.Ux2[ii,jj] = self.P2[ii,jj] / self.Z
            #self.Ux2[ii,jj] = 0

        for ind in index.Ux.minus:
            ii = ind[0] 
            jj = ind[1]
            self.Ux2[ii,jj] = -self.P2[ii+1,jj] / self.Z
            #self.Ux2[ii,jj] = 0

        # update Uy  ===============================
        for ind in index.Uy.field:
            ii = ind[0]
            jj = ind[1]
            self.Uy2[ii,jj] = self.Uy1[ii,jj] - (self._dt/(self._dxyz[1]*self._rho)) * (self.P2[ii,jj+1] - self.P2[ii,jj])
       
        for ind in index.Uy.plus:
            ii = ind[0]
            jj = ind[1] 
            self.Ux2[ii,jj] = self.P2[ii,jj] / self.Z
            #self.Ux2[ii,jj] = 0
        for ind in index.Uy.minus:
            ii = ind[0]
            jj = ind[1] 
            self.Ux2[ii,jj] = -self.P2[ii,jj+1] / self.Z
            #self.Ux2[ii,jj] = 0

        # update P1, Uxy1
        self.P1 = self.P2.copy()
        self.Ux1 = self.Ux2.copy()
        self.Uy1 = self.Uy2.copy()

    def update3(self, index):
        pass

    def plotter(self,s, cb, title, count):
        s.mlab_source.scalars = self.Ux2[1:-1,1:-1]
        time = count/self._fs*1000
        title.text = "sec: %4d [ms]" % time

    def pos2ind(self, position):
        if len(position) == 2:
            return floor((position[0]/self._size.X)*self._num.X), floor((position[1]/self._size.Y)*self._num.Y)
        elif len(position) == 3:
            return floor((position[0]/self._size.X)*self._num.X), floor((position[1]/self._size.Y)*self._num.Y), floor((position[2]/self._size.Z)*self._num.Z)
        else:
            raise ValueError, "Please input correct position." 

# ======================================================
# functions section
# ======================================================

def attributesFromDict(dict):
    """set attributes to this object.
    this code refer to 'PythonCookBook' p284 from O'reilly in Japan.
    """
    self = dict.pop("self")
    for name,value in dict.iteritems():
        setattr(self, name,value)

def wavread(filename, normalize = True):
    """this is like matlab function.
    """
    win = wave.open(filename, 'r')
    buf = win.readframes(win.getnframes())
    x = array(struct.unpack('%dh' % (len(buf) / win.getsampwidth()), buf), float64)
    x = reshape(x, [win.getnframes(), win.getnchannels()])
    if normalize:
        x /= 2 ** (win.getsampwidth() * 8 - 1)
    return x, win.getframerate(), win.getsampwidth() * 8

def wavwrite(x, fs, bits, filename, normalize = True):
    """this is like matlab function.
    """
    win = wave.open(filename, 'w')
    win.setnchannels(int16(size(x, 1)))
    win.setsampwidth(int16(bits / 8))
    win.setframerate(fs)
    if normalize:
        x *= 2 ** (bits - 1)
    buf = array(x, int16).tostring()
    win.writeframes(buf)
    win.close()

def ind2sub(matSize, index):
    """change index number to sub-index.
    input >>>
    matSize: matrix size. [xsize, ysize]
    index: 
    """
    sub=[]
    for ind in index:
        sub.append([ind/matSize[1], ind%matSize[1]])
    return array(sub)

def inpoly(points, nodes, edge = [], TOL = 1.0e-12):
    """inpoly : Point-in-polygon testing
    Determine whether a series of points lie in the bound of polygon
    
    input >>>
    points :
    nodes : 
    edge : 
    TOL : 
    
    output >>>
    cn : 
    on : 
    
    This function based on below's on matlab mFile.    
        Darren Engwirda: 2005-2007
        Email          : d_engwirda@hotmail.com
        Last updated   : 23/11/2007 with MATLAB 7.0
    """
    
    ## error check ##########################################################################
    
    p = points.copy()
    node = nodes.copy()

    try:
        if not isinstance(p, ndarray):
            raise ValueError('P must be an ndarray object.') 
        
        if not isinstance(node, ndarray):
            raise ValueError('node must be an ndarray object.') 
        
        nnode = size(node,0)
                
        if edge == []:   # Build edge if not passed
            temp1 = range(0,nnode)
            temp2 = range(1, nnode)
            temp2.append(0)
            edge = array([temp1, temp2]).T
            
        if size(p,1) != 2:
            raise ValueError('P must be an Nx2 array.') 
            
        if size(node,1) != 2:
            raise ValueError('NODE must be an Mx2 array.') 
            
        if size(edge,1) != 2:
            raise ValueError('EDGE must be an Mx2 array.') 
            
        if (edge.max() > nnode-1) | (any(edge < 0)):
            raise ValueError('Invalid EDGE.') 
            
    except ValueError:
        print 'Please Check the Value'
        raise
    
    ## PRE-PROCESSING ###########################################################################
    n  = size(p,0) 

    # Choose the direction with the biggest range as the "y-coordinate" for the
    # test. This should ensure that the sorting is done along the best
    # direction for long and skinny problems wrt either the x or y axes.
    
    dxy = [max(p[:,0])-min(p[:,0]), max(p[:,1])-min(p[:,1])]  
        
    if dxy[0]>dxy[1]:
        # Flip coordinates if x range is bigger
        p = hstack([p[:, 1:2], p[:, 0:1]])  # swap X- Y 
        node = hstack([node[:, 1:2], node[:, 0:1]]) # swap X- Y 
    
    tol = TOL*min(dxy) 

    # Sort test points by y-value
    ind = argsort(p[:,1])
    x,y = p[ind][:,0:1], p[ind][:,1:2]

    ## MAIN LOOP ###########################################################################
    cn = zeros([n,1], int)    
    # Because we're dealing with mod(cn,2) we don't have
    # to actually increment the crossing number, we can
    # just flip a logical at each intersection (faster!)
    on = zeros([n, 1], int)
    
    for n1, n2 in edge:         # Loop through edges
        # Endpoints     - sorted so that [x1,y1] & [x2,y2] has y1<=y2
        #                       - also get xmin = min(x1,x2), xmax = max(x1,x2)
        y1 = node[n1,1] 
        y2 = node[n2,1]
        
        if y1<y2:
            x1 = node[n1,0]
            x2 = node[n2,0] 
        else:
            yt = y1 
            y1 = y2 
            y2 = yt 
            x1 = node[n2,0] 
            x2 = node[n1,0]
        
        if x1>x2:
            xmin = x2
            xmax = x1
        else:
            xmin = x1 
            xmax = x2 
            
        # Binary search to find first point with y<=y1 for current edge
        if y[0] >= y1:
            start = 0
        elif y[n-1] < y1:
            start = n       
        else:
            lower = 0
            upper = n -1
            for j in range(n):
                start = int(round(0.5*(lower+upper)))
                if y[start] < y1:
                    lower = start
                elif y[start-1] < y1:
                    break
                else:
                    upper = start
        
        # Loop through points        
        for j in range(start, n):
            # Check the bounding-box for the edge before doing the intersection
            # test. Take shortcuts wherever possible!
            Y = y[j]    # Do the array look-up once & make a temp scalar
            if Y <= y2:
                X = x[j]    # Do the array look-up once & make a temp scalar
                if X >= xmin:
                    if X <= xmax:
                        # Check if we're "on" the edge
                        on[j] = on[j] | (abs((y2-Y)*(x1-X)-(y1-Y)*(x2-X)) < tol) 
                        # Do the actual intersection test
                        if (Y<y2) & ((y2-y1)*(X-x1)<(Y-y1)*(x2-x1)):
                            cn[j] = not cn[j] 
                elif Y<y2:   # Deal with points exactly at vertices
                    # Has to cross edge
                    cn[j] = not cn[j] 
            else:
                # Due to the sorting, no points with >y
                # value need to be checked
                break

    # Re-index to undo the sorting
    result = zeros([n,1])
    result[ind] = on
    onIndex = find(result == 1)    
    
    result[ind] = result[ind] + cn
    inIndex = find(result == 1)
    outIndex = find(result == 0)
    
    return inIndex, onIndex, outIndex

    
# [E.O.F]